data {
  int<lower=0> N;                       // Number of data points (observations)
  int<lower=0> N_municipalities;        // Number of municipalities
  int<lower=0> NbDeath[N];              // Observed number of deaths for each data point
  int<lower=0> Population[N];           // Population size for each data point
  int<lower=0> Municipality_code[N];    // Municipality index for each data point (1-based)
  int<lower=0> covid[N];                // Indicator for COVID period (0 = pre-COVID, 1 = COVID)
}

parameters {
  real log_lambda0[N_municipalities];   // Baseline log mortality rate for each municipality
  real log_lambda_c[N_municipalities];  // Log excess mortality rate during COVID for each municipality
  real mu0;                             // Hyperprior mean for baseline log mortality
  real<lower=0> sigma0;                 // Hyperprior standard deviation for baseline log mortality
}

model {
  real lambda_eff[N];                   // Effective expected deaths for each data point

  // Priors for hyperparameters
  mu0 ~ normal(0, 10);                  // Prior for mean of baseline log mortality
  sigma0 ~ normal(0, 10);               // Prior for SD of baseline log mortality

  // Priors for municipality-level effects
  log_lambda0 ~ normal(mu0, sigma0);    // Baseline log mortality per municipality
  log_lambda_c ~ normal(0, 10);         // COVID log excess mortality per municipality

  // Likelihood: Poisson model for observed deaths
  for (i in 1:N) {
    // Compute expected deaths for each data point
    lambda_eff[i] = Population[i] * exp(
      log_lambda0[Municipality_code[i]] + 
      covid[i] * log_lambda_c[Municipality_code[i]]
    );
  }

  NbDeath ~ poisson(lambda_eff);         // Observed deaths follow Poisson distribution
}
